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I. INTRODUCTION 


When formulating and solving an identification problem, it is important 
to have the purpose of the identification clearly in mind. In control 
problems such as control of large space structures, the final goal is to 
design control strategies for a particular structural system. On the other 
hand, there are situation where the primary interest is to analyze the 
properties of a dynamic system, such as stiffness, damping, frequencies, 
etc. The control problem might require a fairly accurate model of the 
system dynamics which will adequately describe the system's motion. 
Fundamentally, one seeks to find a set of parameters that builds a 
mathematical model to best reproduce, according to some criteria, the test 
data. The mathematical model for a linear finite-dimensional dynamic system 
typically includes a state matrix, an output influence matrix and an input 
influence matrix. For flexible structures, damping and frequencies 
constitute the state matrix, mode shapes produce the output influence matrix 
and modal participation factors yield the input influence matrix. In the 
controls field, the process of constructing a model (state space 
representation) from experimental data is called system realization (Refs. 
1-6). The choice of model structure is one of the basic ingredients in the 
formulation of the identification problem. The choice will influence the 
character of the realization problem, the computational effort, the 
possibility for having a minimum order model, etc. The accuracy of 
identified modal parameters from the system realization are affected by the 
choice of the model. 

The purpose of this report is to present methods using experimental 
data to estimate dynamic properties such as damping, frequencies, mode 
shapes and modal participation factors which are referred to as modal 
parameters. The identification of modal characteristics is perhaps the most 
rapidly developing technology in the area of structural system 
identification today. This results from advances in digital data 
acquisition and processing, and the availability of small portable 
minicomputers and microprocessors. Modal surveys which once took months now 
have the potential of being completed in days. The task of modal parameter 
identification is treated in several ways by different researchers (Refs. 7 ~ 
8). There remains considerable controversy about the best technique to use. 
New methods are suggested en masse and, as a result, the field appears to 
look more like a bag of tricks than a unified subject. Therefore, a 
difficult assignment is to organize the bewilderingly large number of ideas. 
Many so-called different methods are in fact quite similar in the sense that 
they are mathematically equivalent but are presented in different ways. 
Thus a unified mathematical framework to treat modal parameter 
identification is needed to achieve some unification of the field. Although 
it is very difficult to get a complete overview of the field which is in 
rapid development, the reader is presented with a basic mathematical 
foundation which provides insight into the field. The objective of this 
report is to offer a unified approach using system realization theory to 
correlate several modal parameter identification methods, particulary for 
flexible structures. 

The tuned sine dwell approach is the earliest method of modal parameter 
identification. The method requires the employment of sine-wave inputs 
having variable frequencies for frequency reponse identification (Refs. 9- 
11), which makes it possible to determine frequencies and, for lightly 
damped systems, damping coefficients. On the other hand, the recently 
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developed methods for modal parameter identification are based on parametric 
models in terms of state equations. Fundamentally, two approaches can be 
classified with respect to the basic elements for construction of the data 
matrix such as the Hankel matrix (Ref. 1). The first approach, which uses 
the impulse response function to construct the data matrix to realize a 
model for modal parameter identification, is called the time-domain 
identification methods (Refs. 12-24). The second approach, which use3 the 
transfer function matrix to realize a model and then identify the modal 
parameters, are referred to as the frequency-domain methods (Refs. 26-29). 
The first part of this report presents the time-domain identification 
methods using system realization theory. The frequency-domain 
identification methods will then be addressed. The relation among several 
existing methods is established and discussed. 


2 



II. TIME-DOMAIN MODAL PARAMETER IDENTIFICATION (Refs. 12-24) 

Many different time-domain methods and techniques in structures field 
were developed, analyzed and tested for modal parameter identification 
(Refs. 12-24). The question arises whether there exists relationship among 
these methods. The answer is positive. Indeed, a unified mathematical 
framework can be developed to present and discuss these methods. The time- 
domain methods for modal parameter identification in the structures field 
can start with the transfer function matrix, which yields Markov parameters 
(generalized impulse response samples). The knowledge of Markov parameters 
makes it possible to construct a Hankel matrix (Refs. 1-6) as the basis for 
the realization of a state space discrete-time model. This section will 
start with the derivation of discrete-time models from the continuous-time 
models which are usually used by structural engineers and follow by the 
basic concept of minimum realization which was developed by Ho and Kalman 
(see Ref. 1). Since the Eigensystem Realization Algorithm (ERA) (Refs. 12- 
15) for modal parameter identification was developed using minimum 
realization theory, it will be presented and discussed first. The 
Polyreference Technique (Refs. 16-18) and the Least Squares Regression 
method will then be derived using the mathematical framework developed in 
deriving the ERA. 


II— 1 . State Equations: Continuous-time and Discrete-time models 
(Refs. 2-6) 

The equations of motion for a finite-dimensional linear dynamical 
system are a set of n 2 second-order differential equations, where n 2 is 

the number of independent coordinates. Let M, D, and K be the mass, 
damping and stiffness matrices, respectively.' The state equations can be 
expressed in matrix notation as 

Mw + Dw + Kw = f(w,t) (1) 


where w , w and w are vectors of generalized acceleration, velocity, and 
displacement respectively and f(w,t) is the forcing function over the period 
of interest at certain specific locations. Eq. (1) can be rewritten as a 
first order system of differential equations in a number of ways. Certain 
reformulations are more suitable for computation than others. One 
reformulation begins with the following definition 

f(w,t) = B f u(t) (2) 


where B f is an n (n = 2n z ) by m input influence matrix. The integer m 

is the number of inputs. Equation (1) can thus be written in a more compact 
form as 
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Ax + Bu 


(3) 



If the dynamic system is excited and measured by the l output quantities in 
the output vector y(t) using sensors such as strain gages, accelerometers, 
etc., a matrix output equation can be formulated as 

y = Cx (4) 

where C is thus an & x n output influence matrix. 

Eqs. (3) and (4) constitute a continuous-time model for finite- 
dimensional dynamic systems. The matrix A in Eq. (3) is a representation 
of mass, stiffness, and damping properties. The input matrix B 
characterizes the locations and type of input u, whereas the output matrix 
C describes the relationship between the state vector x and the output 
measurement vector y. 

Using the initial conditions x(t ) at time t = t gives the solution to 

o o 

the continuous-time equation with an input u(t) 

x(t) = e A ^ fc t o^x(t Q ) + / e A ^ fc 1 ^ Bu ( t ) d-t (5) 

o 

for t 2 t Q . This equation describes the change of state variable x with 

respect to the initial conditions x(t Q ) and the input u(t). It will be 

shown in the following that the evaluation of x(t) at equally spaced 
intervals of time t can be obtained by a discrete-time representation of 
Eq. (5). 

Let the equally spaced times be given by 0, At, 2At,... , kAt,... where 

At is a constant interval. Substitution of t = (k + 1)At and t = kAt into 

o 

Eq. (5) yields 

x[(k+1 ) At] = e AAt x (kAt) + l^ t 1)At e AC(k+1)At " T] Bu( T )d-: (6) 

If u(t) is assumed to be constant over the interval kAt < x < (k + 1)At and 
has the value u(kAt), Eq. (6) with a constant matrix B becomes 

x[(k + 1 ) At] = e AAt x (kAt) + [ {* e^’dx'B] u(kAt) (7) 

where the variable t in Eq. (6) has been changed by letting t’ =(k+1 ) At-x . 
Now, define 

A = e AAt , B = J e AT 'dT'B, x(k+1) = x[(k+1)At], u(k) = u(kAt), (8) 

Note that, for impulse responses, B = B. Eq. (7) can be then expressed in a 
compact form 

x(k+1) = A x(k) + Bu(k) ; k = 0, 1, 2, ... (9) 

and Eq. (4) becomes 
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y(k) = C x(k) 


( 10 ) 


Combination of Eqs. (9) and (10) is the discrete-time representation of a 
dynamical system. This set of equations are the basic formulations for 
system identification of linear, time-invariant dynamical systems, because 
experimental data are discrete in nature. If experimental data are recorded 
in a digital computer, continuous physical measurements will be converted 
into numbers in some format. 

What are the response characteristics of the discrete model, i.e., Eqs. 
(9) and (10)? To observe the response to an impulse in one of the input 
variables, u^O) =1 (i - 1, .... m) and u i (k)=0 (k = 1, 2, ...) are 

substituted into Eq. (9). When the substitution is performed for each input 
element, the results can be assembled into an impulse response function 
matrix Y with dimensions l by m as follows: 

Y(0)=CB, Y( 1 ) =CAB, Y(2)=CA 2 B, ... , Y(k)=CA k B, Y(k+1 ) =CA k+1 B, ... (11) 

This sequence of constant matrices are known as Markov parameters which can 
be obtained from the experimental data through the transfer function or 
impulse responses. The Markov parameters can thus be used as the basis for 
building mathematical models for dynamical systems. 


II-2. Basic Concepts of Kealization (Refs. 1-6) 

The triple of constant matrices [A, B, C], which represents the system 
characteristics can be used to determine the system's response at any of the 
l output points to any input at any of the input points. Such a 
representation is called a realization of the system. Any system has an 
infinite number of realizations which will predict the identical response 
for any particular input. 

Let a new vector z be defined such that 

x = Tz (12) 

where T is any nonsingular square matrix. Substitution of Eq. (12) into 
Eqs. (9) and (10) yields 

z(k+1 ) = T -1 A T z(k) + T _1 Bu(k) ; k - 0, 1, 2, ... (13) 

y(k) = CTz(k) (14) 

It is obvious that the effect of input u(k) on y(k) will be the 3ame for 

this new system of equations (13) and (14). Thus, the triple [T ^ AT , T ^ B , 
CT] will also be a realization for the same system. The predicted responses 

using the realization [T 1 AT, T 1 B, CT] will be identical to those predicted 
using [A, B, Cj. Because there exists an infinite number of nonsingular 
matrices T, there are an infinite number of such realizations. 

By minimum realization is meant a model with the smallest state space 
dimension among systems realized that has the same input-output relations. 
All minimum realizations have the same set of eigenvalues which are 
parameters of the system itself. 
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Assume that the state matrix A of order n has a complete set of 

o 

linearly independent eigenvectors ( 4 ^ , 4 > 2 ,..., 4 / ) with corresponding 

o 

eigenvalues (X^ , X 2 ) which are not necessarilly distinct. Define 

0 

A as the diagonal matrix of eigenvalues and 41 the matrix of eigenvectors, 
i .e. 


A = diag. ( X.j , X 2 X n ) (15) 

o 

and 

4) = Lvl> 1 , \l> 2 4> n 3 (16) 

o 

The triple [A, B, C] can then be transformed to the realization [A, 4 > ^ B , 
C 4 j ] by choosing T = ip. The diagonal matrix A contains the information of 

modal damping rates and damped natural frequencies. The matrix i|> ^B is 
called initial modal amplitudes and the matrix C 41 mode shapes. All the 
modal parameters of a dynamic system can thus be identified by the triple 

[A, 4 > ^ B , C 4 >]. The desired modal damping rates and damped natural 

frequencies are simply the real and imaginary parts of the eigenvalues A, 
after transformation from the discrete-time domain to the continuous-time 

domain using the relation A = ln(A)/At. 


Time-domain analysis for identification of modal parameters begins by 
forming the generalized (r+1) by (s+1) Hankel matrix, composed of the Markov 
parameters (see Eq .11). 


H(k) 


Y(k) Y(k+1 ) ... Y(k+s) 

Y(k+1 ) Y(k+2) ... Y(k+s+1 ) 


Y(k+r) Y(k+r+1 ) ... Y(k+r+s) 


(17) 


If r+1 2 n and s+1 2 n (the order of the system), the matrix H(k) is of 
0.0 

rank n Q . To confirm this point, observe from Eq. (11) that 


H(k) = V A k W (18) 

r s 

where 

W = [ B AB A 2 B ... A S B ] (19) 

s 

is called the observability matrix, whereas the block 

matrix W is called the controllability matrix. If the order of the 
. s 

system is n Q , then the minimum dimension of the state matrix is n Q x n Q . 


C 

CA 


CA 


anc 


CA 


The block matrix V 
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If the system is controllable and observable, the block matrices V and W 

r s 

are of rank n . Therefore the Hankel matrix is of rank n by Eq. (18). 

o o M 

Based on the properties of the Hankel matrix composed of the Markov 

parameters (pulse response function), several methods for modal parameter 

identification are discussed in the following sections. 

II-3. The Eigensystem Realization Algorithm (ERA) (Refs. 12-15) 

The basic development of the time-domain (state-space) concept is 
attributed to Ho and Kalman (Ref. 1) who introduced the important principles 
of minimum realization theory. The Ho-Kalman procedure uses the Hankel 
matrix (Eq. 17) to construct a state-space representation of a linear system 
from noise-free data. The methodology has been recently modified and 
extended to develop the Eigensystem Realization Algorithm (Refs. 12-15) to 
identify modal parameters from noisy measurement data. 

In contrast to classical system realization methods which use the 
Hankel matrix (Eq. 17), the ERA algorithm begins by forming a block data 
matrix which is obtained by deleting some rows and columns of the Hankel 
matrix (Eq. 17), but maintaining the first block matrix intact. Furthermore, 
the standard ordering of entries in the Hankel matrix does not need to be 
maintained. 

T T T T 

Let B = [b, , b^,..., b ] and C = [c. , c„ , . . . c„ ] where the column 

12m 1 2 ft 

vector b^ is the control influence vector for the ith control input and the 

row vector c. is the measurement influence vector for the jth measurement 
J 

sensor. Denote that columns of B r (i = 0, 1, ..., n) and the rows of C. (j 

L J j 

= 0, 1, .... 5) as arbitrary subsets of (b, , b„,..., b ) and (c. , c_,..., 

1 d m i d 

T 

c^), respectively. If the column index 1^ is defined as the set of integers 

(1, 2, ..., m) , the column index I ^ ( i =1,2, ..., n) will be any arbitrary 

subset of Iq. On the other hand, if the row index J Q is the set of 

integers (1, 2, ..., 1) , the row index J.(j = 1, 2,..., £) will be any 

J 

arbitrary subset of J Q . The ERA block data matrix can then be expressed by 


S + + 1 

H(k) = [Y (s +k+t )]; Yj (s.+k+t ) = A i JB 

j i J i J j i 


( 20 ) 


where s Q = t Q = 0, and s^ and t., are arbitrary integers. When i = j = 0, 
Y. r (k) = Y(k) = CA k B. 

-Vo 


The ERA block data matrix (Eq. 20) allows one to include only good or 
strongly measured signals without losing any capability. This is useful 
since some measurement data may be noisier than others or sensors may 
malfunction during the test. The advantage of this capability is the 
potential to minimize the distortion of the identified parameters caused by 
noise. A judicious choice of data and its proper arrangement in the block 
matrix H(k) can also be used to minimize the computational requirements of 
the method. For example, the columns of H(k) may be made as independent as 
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possible by properly selecting the data samples to use as entries of the 
matrix. This effort could substantially reduce the order of the matrix for 
large problems. For sufficiently low noise data, the order can be the same 
as that of the true system state matrix A. This fact results from 
examination of the controllability and observability matrices, to be 
discussed next. 


From Eqs. (20), it can be shown that 

1C 


H(k) = V A W ; V_ * 
5 n 5 


C A S 1 
J 1 


C A S S 
J 5 


; and u - [B A C t B, ... » t nB r J (21) 


where 


V and W 
5 n 


Assume that 


W H # V P = I 
n S r 


are generalized observability and controllability matrices, 
there exists a matrix satisfying the relation 

( 22 ) 


where I is an identity matrix of r. It will be shown that the matrix H # 
r 

plays a major role in deriving the ERA. What is H^? Observe that 

H(0)H # H(0) = V r W H # V r W = VW = H(0) (23) 

5 n £ n 5 n 


The matrix is thus the pseudoinverse of the matrix H(0) in a general 
sense. 


The ERA process starts with the factorization of the block data matrix 
(Eq. 20), for k=0, using singular value decomposition (Ref. 29): 

H( 0) = P D Q T (24) 


where the columns 
rectangular matrix 


of matrices 
D 


D = 


0 


0 

0 


P and Q are 
with = diag 


d^] and monotonically non-increasing d^i-1, 2, 


orthonormal and D 


C d i , d_, ... , 


n’ n+1 


. . . ,r) 


is a 


a d 2 J 


> d > d > 

n n+1 


> d > 0. 
r 


Next, let P and Q be the matrices formed by the first r columns of P and 
r r 

Q, respectively. Hence the matrix H(0) and its pseudoinverse become 


H(0) = PDQ 
r r r 


where P P 
r r 


T 

I n = 

r r r 


( 25 ) 


and 
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( 26 ) 


H 


# 


Q D _1 P 
r r 


T 

r 


Eq. (26) can be readily proved by observing Eq. (23). 

Define 0^ as a null matrix of order SL, I an identity matrix of order 

t, and eJ = [I 0 ...0 ]. Using Eqs. (20), (21), (22), (25), and (26), a 
minimum order realization can be obtained as follows: 


Y(k) 


eT H(k)E = eJv A^W E 
i m 5 n m 

eIv [W H # V ]A k [W hV]W E 
x. S n S n 5 n m 

eTh(o)[q d” 1 pJ]v a^w [q d" 1 P^]H(0)E 
x, r r r t, n^rr m 

e!h( 0) Q D~ 1/2 [D~ 1/2 P T H(1)Q D _1 72 ] k D~ 

Z rr r r rr r 

eTp D 1/2 [D- 1/2 P T H(1)QD- 1/2 ] k D 1/2 Q T E 
£ r r r r rr r rm 


(use Eqs. (20) & (21)) 

(use Eqs. (22) 

(use Eqs. (21 ) & (26) ) 

P^H ( 0 ) E (use Eq. (25)) 
r m 

(27) 


This is the basic formulation of realization for the ERA. The triple 


A = D _1/2 P T H(1)Q d" 1/2 
r r r r 


d 1/2 q t e 

r r m 


C = 


T_ 1 /2 
E.P D 
l r r 


(28) 


is a minimum realization. The order of the matrix A is r which is equal 
to n Q (the order of the system) for sufficiently low noise data. When the 

matrices P and Q are obtained through other factorization methods such 
r r 

that P T P * I and Q T Q * I , Eq. (28) is still valid if the matrices P^ and 
r r r r r r r 

T // # 

are replaced by P^ and respectively. 


Due to measurement noise, nonlinearity, and computer roundoff, the 
block matrix H(k) will usually be of full rank which does not, in general, 
equal the true order of the system under test. It should not be the aim to 
obtain a system realization which exactly reproduces the noisy sequence of 
data. A realization which produces a smoothed version of the sequence, and 
which closely represents the underlying linear dynamics of the system, is 
more desirable. Several accuracy indicators have been investigated for 
quantitatively partitioning the realized model into pure (principal) and 
noise ( perturbational ) portions so that the noise portion can be 
disregarded. Two principal indicators now available are the singular values 
of the block data matrix and a parameter referred to as Modal Amplitude 
Coherence. The number of retained singular values determines the order of 
the realization, and Modal Amplitude Coherence is used to assess the 
resulting degree of modal purity. 

If Eqs. (21) and (22) are examined as a whole the equality 

H( 0) = V W = [P d]/ 2 ] [DyV] (29) 

£ n r r r r 

defines the controllability and observability Grammians as 
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( 30 ) 


T 

W W 
n n 


D and vTv = D . 
r E, 5 r 


The fact that the controllability and observability Grammians are equal and 
diagonal implies that the realized system [A,B,C] is as controllable as it 
is observable. This property is called an internally balanced realization. 
It means that the signal transfer from the input to the state and then from 
the state to the output are similar and balanced. 

Some singular values, say d n+ ^,...d , may be relatively small and 

negligible in the sense that they contain more noise information than system 
information. In other words, the directions determined by the singular 
values d n+ j . .d r , have less significant degrees of controllability and 

observability relative to the noise. It would be unwise to require a 
realization including these directions. The reduced model of order n 
after deleting singular values d d is then considered as the 

robustly controllable and observable part of the realized system. Reference 
13 provide the mathematical framework for establishing the relationship 
between these accuracy indicators and the characteristics of the noise. 


An Alternate ERA Formulation for Implementation in a Small Computer . 

For a large flexible dynamic system under test, the number of sensors 
may be on the order of one hundred which correspondingly can result in the 
number of rows for the block matrix H(k) ( Eq . 20)) being significantly 
larger than the order of system. This then requires a large amount of 
computer storage and computational efforts. To circumvent this difficulty, 
an alternate formulation is derived in the following. 

Let Tp be a rectangular matrix of rank at least r which operates on the 
rows of the block matrix H(k) (Eq. 20)) such that 

H(k) = T H(k) = [T V r ]A k W = V A k W (3D 

r r C n 5 n 

The operator of on H(0) may include: (1) multiplying a row by a nonzero 

number; (2) interchanging two rows and (3) adding the product of one row to 

another row. As shown in Eq. (22), the pseudoinverse H of the matrix H(0) 
satisfies the relation 

W H # V r = I (32) 

n 5 

Now, factorize the data matrix H(0), using singular value decomposition, 

H( 0) = P 5 where pl? n = I - oJq (33) 

nnn nn n nn 

and thus 

H # = QD' 1 ^ (3D 

n n n 
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The integer n is determined by retaining only the non-negligible singular 
values. Following the similar procedure shown in Eq. (27), a minimum order 
realization can be obtained as 


V A k W = VTW H # V r ]A k [W H # V.]W 
5 n 5 n 5 n 5 n 

= H(0)H # [V^.A k W Ti ]H # H(0) 

- [P D^KD^Vv A k W Q d" 1/2 ][5’ /2 Q T ] 

n n n n g n n n n n 

= [p d 1 / 2 ][d” 1 / 2 p t h(i)q d _1/2 ] [d 1 / 2 q t ] 
n n n n.nn n n 


(use Eq. (32)) 

(use Eq. (31 ) ) 

(use Eqs. (33) & (3*0) 

(35) 


The triple 


v = p d 1/2 , a = 5 1/2 p t h(i)q d” 1/2 , W = 5)/ 2 q t 

5nn n n . nn qnn 


(36) 


is a minimum realization of order n in which the. state transition matrix A, 
the controllability matrix and the modified observability matrix after 
the row operation by T^ are realized. 

From the definition of the controllability matrix in Eq. (23) and 
T 

the matrix E = LI 0 ... 0 ], the input matrix B can be expressed as 

m mm m 


B 


5y 2 Q T E 

n n m 


(37) 


which takes the first m columns of the matrix W^. To recover the output 

matrix C, a simple computation can be derived by observing from Eqs. (21) 
and (36) that 

H(0) = V_W = V'D'/ 2 *: (38) 

which produces the original observability matrix as 

V = H(0)Q 5" 1/2 (39) 

5 n n 

The output matrix C can thus be obtained from the observability matrix by 

c = eJv = eJh(o)q D ' 1/2 ( *10) 

IE, l n n 

As a result, the triple 

A = 5 n' /2 ' P 5 (,) Vn 1/2 » B - B y 2 Qn E m * C = E o H(0) Q n 5 n 1/2 (lj1) 

n n.nn nnm x, nn 

is a minimum realization of order n for the system. 

The advantage of this formulation is that the modified matrix H(0) can 
represent a reduced-in-size version of the block data matrix H(0). The 



number of rows for the matrix H(0) can be set to a size only moderately 
larger than the order of the system without loosing the information of 
parameters in the data. For example, if the number of rows for the matrix 

H(0) is set to r, all the measurements represented by a number larger than 
r can be added to a particular row or distributively to the first r rows. 

In doing this row operation, the matrix H(0) will be reduced in size, and 
also the measurement data may be averaged to reduce the noise level of the 
data. Note that all the measurements may be added together and then shifted 

to form the matrix H(0). The number r can not be chosen smaller than twice 
the anticipated number of system modes, or some modes would not be 
identified. If the controllability and observability Grammians are examined 
for this realization, it is seen that the realization is not internally 
balanced. Thus, if the system based on this realization is reduced for 
controller design, the reduced system may not be as robust as the one based 
on the realization of Eq. (28). 


II-4. The Polyreference Technique (Canonical-Form Realization) 

(Refs. 2, 6, 16-18) 

During the past two decades, several algorithms for the construction of 
a canonical-form representation of linear systems have appeared in the 
controls literature (Refs. 2 & 6). Controls researchers are mainly 
concerned with, for an example, determining a passive or an active network 
that has a prescribed impedance or transfer function. Although techniques 
of canonical-form realization are available in the controls literature, 
direct application to modal parameter identification for flexible structures 
has not been addressed. Among several available methods for canonical-form 
realization, one based on Hankel matrices will be addressed in this section. 

Recently in the structures field, a method, similar if not identical, 
to a canonical-form realization, was developed in Refs. 16-18 using 
frequency-response functions for identification of modal parameters from 
multi-input and multi-output measurement data. The method is referred to as 
the Polyreference technique. Mathematical background of the Polyreference 
technique will be presented in this section using a new approach which 
provides insights of correlation between the Polyreference technique and the 
Canonical-Form Realization. 


Form the (r+1)x(s+1) block Hankel matrix 


H(0) 


Y(0) Y(1) 

Y ( 1 ) Y(2) 


Y(s) 
Y(s+1 ) 


(42) 


Y(r) Y(r+1) ... Y(r+s) 


where r and s are integers which are chosen to be larger than the order 
of the system. Using the singular value decomposition, find the nonsingular 


matrices 

P and Q such that 




H(0) 

= pdq t = IP P 1 

n o 

ID 

n 

°'| [Q n V T 

(43) 



L° 

0 J 



1 2 



where D n is a diagonal matrix containing monotonically non-increasing non- 

negligible singular values. The integer n is determined by the 
characteristics of the system noise as discussed in Refs. 13 and 29. All 
singular values numbered after n are considered as zero values. The 
matrices and P q denote respectively the first n and the last 

remaining columns of the orthonormal matrix P. Similarly, and Q q 

denote respectively the first n and the remaining columns of the 
orthonormal matrix Q. 


Now observe, from the definition of Markov parameters, that 


H(0 ) = VW; V = 


C 

CA 

CA 1 


and W = [B AB ... A B] 


(44) 


where V and W are observability and controllability matrices 

respectively. Since matrices 
I, Eqs. (43) and (44) lead to 


T T 

respectively. Since matrices P and Q are orthonormal, i.e., P P = Q Q = 


p t h(o)q = p t vwq = 

~ P l VWQ 
n n 

p t vwq 
n o 


1 

O 

c 

Q 


p t vwq 

o n 

p t vwq 

o 0 


1 

O 

o 


which yields 


P VWQ = D , 
n n n 


P T VWQ = 0, p t vwq = 0 
no on 


and 


P VWQ = 0. 
o o 


(45) 


(46) 

(47) 

(48) 


Note that the ERA algorithm is developed using the matrices P^, Q n and D n as 

shown in Eq. (46). If the system is assumed controllable and observable, 

each of the five matrices P , V, W, Q and D is of rank n. This means 

n n n 

T 

that the ranks of matrices P V and WQ are n. Thus Eqs. (47) and (48) 

n n 

imply 


P T V =0 and WQ =0 (49) 

o o 

The matrix P q provides the left orthonormal basis for the null subspace 
which is orthogonal to the observability matrix, whereas the matrix Q q 

gives the right orthonormal basis for the null subspace which is orthogonal 

to the controllability matrix. Now partition the matrices P and Q as 

oo 
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( 50 ) 


P T = C pT p T 

o oO ol 


P T j and Q T = [Q T n Q 1 , . . . Q T ] 
or o oO ol os J 


Substitution of Eqs. (44) and (50) into Eqs. (49) yields 
r 

* *1* 1 ip T rp *>> 

E P .CA = P . C + P ,CA + ... + P CA = 0 
. . oi oO ol or 

1=0 


E A BQ . = BQ + ABQ + 
i = 0 01 00 0 


+ A BQ =0 
OS 


(51a) 


(51b) 


Eqs. (51) are the basic formulation for the Polyreference technique and 
the Canonical-Form Realization. In fact, Eq. (51a) is the basis for an 
observable canonical-form realization whereas Eq. (51b) is the basis for a 
controllable canonical-form realization (see Ref. 4, pp. 321, Problems 6-19 
and 6-21). Both equations (51) should produce the same results. The 
question arises as to whether Eq. (51a) is more favorable than Eq. (51b) or 
vice versa. The answer is given in the following. 


Observe that each submatrix P (i=0,...,r) must have more columns than 

rows (the number of outputs i ) . Similarly, each submatrix Q Qi (i = 0, ..., 

s) has more columns than rows (the number of inputs m) . Form square 

matrices P . of order l and Q . of order m respectively from the last 
oi oi 

% columns of the matrix P . and the last m columns of the matrix Q . , and 

oi oi 

rewrite Eqs. (51) such that 


and 


r [p T ]~ 1 p T c a 1 = c A r = r p T c a 1 

1-0 ° r 01 i=0 01 


- V a 1 B Q oi [Q os ] _1 = A S B = A V A 1 B Q o . 

i=0 1=0 


with 


- T -T-1-T - - - -1 

P = [P ] P ! and Q . = Q . [Q ] . 
oi or oi oi oi os 


(52a) 

(52b) 


Eqs. (52) can be rearranged into companion matrix form as 


c 


0 0 ... 0 


" C 

CA 


0 0 I ... 0 


CA 

: 

A = 

J J J l 


J 

CA r ' 2 


0 0 0 ... I 


CA r " 2 

_ r-1 


~T -T -T -T 


r-1 

CA 


_ _P o0 " P o1 " P o2 P o(r-1)- 


_ CA . 


(53a) 
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and 


A [B AB ... A S_1 B ] = [ B AB ... A S_1 B ] 


0 0 


I 0 
m 


. 0 Q 


oO 


0 Q 


ol 


0 I ... 0 Q„ 
m o2 


0 0 ... I Q . , . 

m o(s-1 ) 


(53b) 


The matrix I. is an identity matrix of order J, and the matrix I is an 
x. m 

identity matrix of order m. Now it is claimed from Eq. (53a) that the 

triple 



1 

o 

o 

CH 

0 

1 


" Y(0) “] 



o 

• 

• 

• 

o* 

O 

o 


Y ( 1 ) 


A = 

0 0 0 ... I 

-T -T -T ~T 

, B = 

Y(r-2) | 

• 

, c = [i a 0 ... 0 0 ] 


-P .p‘ -p 1 -p‘ 

oO ol o2 "• o(r-1 ) J 


[_Y(r-1)J 

(54a) 


is an £,r-dimensional realization of the system, indeed, it can be readily 
verified that 

Y(0) = CB, Y( 1 ) = CAB, .... Y(r) = CA r §. 


Because of the structure of A and C, it is easy to show that the 
realization is observable. However, it is not in general controllable. 
This realization is called an observable canonical-form realization. It is 
not a minimum realization, because, in general, it is not both controllable 
and observable. 


Similarly, from Eq. (53b), it can be verified that the triple 



0 

0 

... 0 

^oO 


r i i 

m 




I 

m 

0 

... 0 

«01 


0 



A = 

0 

I 

m 

... 0 

^o2 

, B = 

0 

, C = LY(0) Y ( 1 ) 

... Y ( s-2) Y(3-1)] 


0 

0 

... I 

m 

Q o(s-l) 


_ 0 


(54b) 


is an ms-dimensional realization of the system. This is a controllable 
canonical-form realization which is not in general observable. Again this 
realization is not of minimum order. 

Eqs. (5*0 can be reduced to a minimum order realization by applying the 
reduction procedure shown in Ref. 6, Chapter 5. However, the canonical form 
will be destroyed after the application of the reduction procedure. 
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The order of either observable or controllable canonical-form 
realization, i.e., £r or ms, is required to be equal to or larger than the 
order of the system. From a computational point of view, one should choose 
the one with smaller dimension to work for modal parameter identification. 
The numerical problem for the ei gensol ut i on of the canonical-form 
realization can be solved in various ways. A technique suitable and 
efficient for mini-computer systems has been implemented and shown in Ref. 
31 . It should be remarked that only a subset of eigenvalues in the realized 

state matrix A (see Eq. (54)) belongs to the actual state matrix A, since 

the matrix A is generally oversized for multi-input and multi-output 
measurements. 

The method used here to obtain the canonical-form realization is 

different from that shown in Ref. 16. Orthonormal matrices P and Q are 

o o 

computed through the application of the singular value decomposition to 
realize a companion-form state matrix. Since orthonormal matrices are very 
close to identity matrices, this thus generates a computationally well- 
behaved canonical-form realization. 


II-5. An Alternate Method for the ERA and the Polyreference Techniques 
(Ref. 19) 

A minimum order of canonical- form realization is generally impossible 
for multi-input, multi-output systems due to the constraint that the 
realized state matrix is a companion form. If the constraint is released, a 
minimum order realization can be obtained from Eq. (53). 

In view of Eqs. (43) and (44), the controllability and observability 
matrices can be expressed by the following equations 

1 /? 1 /2 T 

V = P D and W = D Q (55) 

n n n n 


Define 0 as a null matrix, I. an identity matrix of order ft, r and E. 

x,r x,r 

[I. 0] of dimension Ir x X,(r+1). 

x.r 


Hence, Eq. (53a) can be written with the aid of Eq.(54a) as 


[ VV D r A ■ j£E tr P n ]D n /2 ’> A ' %' / 2 [E *r P / A[ VrX'" < 56 > 


where # means pseudoinverse. This is a minimum realization of order n. To 
compute Eq. (56), a simple procedure can be developed as follows. Define 

as a shift operator which shifts l columns of a matrix, for an example, 

o„E„ = [0 I„ ]. In view of the definition of the shift operator o„, the 
£ £r xr x, 

companion matrix A in Eq. (54a) and Eq. (55), Eq. (56) becomes 


A 


D 1 7 ^CE. P l^Co.E P jD^ 7 ^ 
n x,r n l x,r n n 


(57) 


Here, [o E P ] simply means the matrix obtained by deleting the last l 
l ir n 

rows of the matrix P and [E. P ] represents the matrix obtained by 

n 36 n 
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deleting the first H rows of the matrix P . This equation was first 

presented in Ref. 19. Since P n is an orthonormal matrix, a special and 

efficient procedure can be developed to compute the pseudoinverse of the 
matrix E^P using the matrix inversion lemma (Ref. 19). 


Similarly, an equation for the state transition matrix A can be derived 
from the controllable form (Eqs. (53b)) as 


„ J/2 r ^ ,T r 0 ,T// -1/2 

A = D [o E Q ] [E Q ] D 
n m ms n ms n n 


(58) 


where [E Q ] means the matrix obtained by deleting the first m rows of the 
ms n 

matrix Q and [o E Q ] represents the matrix by deleting the last rows of 
n m ms n 

the matrix Q . 

n 

Eqs. (57) and (58) can also be derived by the ERA procedure. Let an 
oversized Hankel matrix H be formed such that 


H( 0) = E^ H and H(1) = o^E^H (59) 

The Hankel matrix H(0) is formed by deleting the last H rows of the Hankel 
matrix H whereas the Hankel matrix H(1) is obtained by deleting the first H 
rows of the matrix H. 

Find the orthonormal matrices P and Q , and a diagonal matrix D such 

n n n 

that 

H = P D (60) 

n n n 

Eq. (59) thus becomes 

H( 0) - [EPlDQ^ and H(1) = Co.E, P ]D (61) 

Hr n n n H Hr n n n 

Substituting this equation into the ERA basic formulation (^1) and noting 
T 

that Q Q = I yields the triple 
n n n J 


A = D _1/2 [E P ] # [a. E P ]D 1/2 , B = D 1/2 Q T E , C - eTp D^ 72 (62) 

n' Hr n H Hr n n n n m H n n 

The state transition matrix A is indeed identical to that in Eq.(57). 

Similarly, if 


H( 0) = HE T = P D [E Q ] T and H(1) 
ms n n ms n 


H[o E ] T = PD[oE a Qj T (63) 
m ms n n m ms n 


where the Hankel matrix H(0) is obtained by deleting the last m columns of 
the oversized Hankel matrix H and the Hankel matrix H(1) is obtained by 
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deleting the first m columns of the same matrix H. Substitution of Eq. ( 63 ) 

T 

into the ERA basic formulation, with the aid of P P = I , produces the 

n n n 

triple 


d' /2 U E « ] T [E « ] T V ’ /2 , B= D V2 Q T E 
n m ms n ms n n n n m 


T 1 /2 

c " Wn (64) 


The state transition matrix A is again identical to that in Eq. (58). 

Realizations (62) and (64) preserve the same features as for the ERA, 
including a good numerical performance, internal balancedness, and 
flexibility in determining order-error tradeoff. Based on formulations (62) 
and (64), a close link between the ERA and the Polyreference techniques is 
established through the singular value decomposition and the generalized 
Hankel matrix. 


II-6. Least Squares Regression Techniques (Refs. 20-24) 

The least square regression technique for a discrete-time dynamic model 
has been derived and used for system identification for more than two 
decades (see Ref. 20, Chapter 5, pp. 97-99). The same technique was 
rederived and further developed for modal parameter identification in the 
structures field (see Ref. 23-24). Here, the least squares regression 
technique will be formulated using system realization theory which provides 
a good basis for the comparison with other methods. 

In view of Eqs. (21) and (22), the measurement function Y(k) can be 
obtained through either of two other algorithms as follows: 

Y(k) = E^H(k) E = e”Jv _A k [W H # V_] W E m 
l m IE. n 5 

- e„[v a w H # ] k v c w E 
i i n 5 n m 

= E^[H(1)H # ] k H(0)E m (65) 

i m 

or 

Y(k) = E^H(k)E = E^V.CW H # V.]A k W E m 
l m x, t, n t, pm 

= E^vw [H^V A W ] k E 

1 ( n 5 n m 

= e]h( 0)[H # H(1 ) ] k E (66) 

x, m 

Hence, the triple [H(1)H # , H(0)E , E,] or the triple [H # H(1), E , eJh(O)] is 

m X# m x< 

a realization. The matrix H(1)H # or H^H(1) constitutes the basis for the 
least square regression technique (see Ref. 20, pp. 97-99). 

The matrix H # is the pseudoinverse of the matrix H(0). For a single 
input, there exists a case where the rank of H(0) equals the column number 
of H( 0) , then 

# T -1 T 

H = [H( 0) n( 0) ] H(O) 1 
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On the other hand, there exist a case for a single output where the rank 
equals the row number, then 

H # = H(0) T [H(0)H(0)V 

The matrix H(1)H^ has been used in the structural dynamics field to identify 
system modes and frequencies (see Refs. 23-24). This is a special case 
representing a single input which cannot realize a system that has repeated 
eigenvalues, or using sufficiently low noise data unless the system order is 
known a priori. 

These realizations (Eqs.(65) and (66)) are not of minimum order, since 
the dimension of the state vector, x, is the number of either columns or 
rows of the matrix H(0) which is larger than the order of the state matrix A 
for multi-input and multi-output cases. Examination of Eqs. (65) and (66) 
reveals that these two equations are special cases of ERA. Eq.(65) is 
formulated by inserting the identity matrix (Eq. (22)) on the right-hand 
side of the state transition matrix A. On the other hand, Eq. (66) is 
obtained by inserting the identity matrix (Eq. (22)) on the left-hand side 
of the state transition matrix A. However, the ERA is formed by inserting 
the identity matrix (Eq. (22)) on both sides of the state transition matrix 
A as shown in Eq. (27). Because of these differences, the least squares 
regression method does not minimize the order of the system. 
Mathematically, if the singular value decomposition technique or other rank 
detection technique is not included in the computational procedures, the 
realized triple obtained from Eq. (65) or (66) cannot be numerically 
implemented, unless a certain degree of artificial noise and/or system noise 
is present. Noise tends to make up the rank deficiency of the data matrix 
H ( 0 ) . Since noise contamination is not quaranteed, the least squares 
regression technique (Eq. (65) or (66)) is subject to uncertain results. 
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III. FREQUENCY-DOMAIN MODAL PARAMETER IDENTIFICATION 

Time and frequency are two fundamental bases of description for linear 
dynamic systems (Ref. 25). For example, given a single input and single 
output linear dynamic system, a three-dimensional space can be constructed 
for the output response with amplitude as one axis, and time and frequency 
as the other axes. A sinusoidal time history for each individual frequency 
(mode) can be treated as a projection on the time plane, existing at some 
distance from the origin. This distance is measured along the frequency 
axis. Similarly, the dynamic output response has a projection onto the 
frequency plane. This projection takes the form of an impulse with an 
amplitude. The position of the impulse coincides with the corresponding 
frequency. Summing multiple time plane projections produces the time 
history of the dynamic response. Similarly, connecting all the frequency 
components in the frequency plane yields the spectrum. The duality of the 
time and frequency description of the dynamic response for a linear system 
becomes evident. 

This section exposes the close conceptual connection between time- 
domain and frequency-domain approaches to identification of modal parameters 
for linear dynamic systems. To identify modal parameters including damping 
ratios, frequencies, mode shapes, and modal participation factors, many 
methods have been developed using classical frequency spectra (transfer 
function) analysis or sophisticated time-domain methods. For many years, 
the multi-shaker sine-dwell approach was the primary experimental method 
employed in aerospace structural mode surveys. The major advantage is that 
each mode is experimentally isolated by careful tuning of a series of 
shakers and verified as part of the experimental process. The disadvantage 
lies in the long elapsed testing time and the requirement for considerable 
test engineering expertise. 

The current section presents some recently developed frequency-domain 
techniques. All of the methods will be correlated using system realization 
theory. There are a number of important features in the frequency-domain 
analysis, including, overlap averaging and zooming (Ref. 25). The overlap 
averaging is used to smooth the transfer function, while the zooming is used 
to concentrate all the spectral lines into a narrow band in the frequency 
range of interest. 

The transfer functions are basic elements for frequency-domain 
techniques. This section starts with the discussion of the transfer 
function characteristics. It then follows by the presentation of 
Eigensystem Realization Algorithm in frequency domain and the Polyreference 
technique in frequency domain, for the purpose of comparison with time- 
domain identification techniques. The conceptual connections between time- 
domain and frequency-domain approaches are exposed and discussed. 


I I I - 1 . Characteristics of the Transfer Functions. 

As pointed out in the last section, the linear, constant, finite- 
dimensional dynamical systems can be represented by the continuous-time 
model 


x = Ax + Bu (67) 

y - Cx (68) 
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Taking the Laplace transform of Eqs. (67) and (68) yields 

y(s) = C[sl - A] 1 (x(t Q ) + Bu(s) ) (69) 

where I is an identity matrix and s is the parameter of the Laplace 
transform. For x(t ) =0, 

y(s) = C[sl - A] _1 Bu(s) (70) 

This equation shows the direct relation between the input spectrum u(s) and 
the output spectrum y(s). The function defined by 

Y(s) = C[sl - A] _1 B (71) 

is called the transfer function matrix of dimension Z (the number of 
outputs) by m (the number of inputs). 

Let Y(t) be the impulse response function matrix 

Y(t) = C e At B (72) 

which is then related to the transfer function matrix Y(s) by 

f OO r 00 — 

Y(s) = L[Y(t) ] = j Q Y(t)e~ St dt = J 0 Ce At § e~ St dt = C[sl - A]" 1 B (73) 
or 

Y(t) = l"![Y(s)] = j Y(s)e st ds = J C[sl - A]" 1 B e St ds = C e At B (74) 

where the eigenvalues of A are assumed to have negative real parts and L ^ [ 
] is the inverse of the Laplace transform operator L[ ]. The Laplace 
transform of the ith derivative of the impulse response function matrix 
becomes 

/•OO - , CO - 

L[ Y 1 ( t) ] = J 0 CA 1 e At B e~ 3t dt = CA 1 J Q e At B e" St dt = C^EsI-A]" 1 ^ = CA 1 G(s)B 

= s X Y(s) - s 1_1 Y(t=0) - s 1_2 Y (1) (t=0)- ... - Y (l_1) (t=0) 

= 3^(3) - l Y (k) (t=0) s i " 1 “ k ; i = 1. 2, ... (75) 

k=0 

where Y^^= d i Y(t)/dt i and G(s) = [sI-A] 1 

Let the equally spaced time be given by 0, At, 2At, ..., kAt, ... where 
At is a constant interval. The numerical Fourier transform of the impulse 
response function Y(jw) has the approximation as follows: 
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Y(joj) = j Y(t)e~ ju>t dt - l Y(kAt)e" Ju)kAt At = l Ce AkAt Be _J “ kAt At 

k=0 k=0 


= l CA Be 
k=0 


k -jmkAt 


At = I Y(k)e 
k=0 • 


-jwkAt 


At = CG(ju)B 


( 76 ) 


where A=e AA ^(Eq. 8), G(ja)) = £ A k e '^ a)kA ''At and Y(k) = CA k B (Eq. 11). The 

k=0 

Fourier transform of the ith derivative of the impulse response matrix can 
correspondingly be approximated by 


(Jw^YCju)- l Y (k) (t=0) 
k=0 


(jm) 


i— 1 — k 


l CA 1 e AkAt Be Ju3kAt At=CA 1 G(jio)B 
k=0 


(77) 


Eqs. (76) and (77) constitute the bases for the frequency-domain parameter 
identification methods which will be presented in the following. 


III-2. The Eigensystem Realization Algorithm in Frequency Domain (ERA-FD) 
(Ref. 26) 

This subsection presents a brief review of a recently developed 
technique for modal parameter identification which is referred to as the 
Eigensystem Realization Algorithm in Frequency-Domain (ERA-FD). Similar to 
the Eigensystem Realization Algorithm in time-domain, the ERA-FD starts with 
a complex data matrix formed by transfer functions and corresponding shifted 
transfer functions. Using singular value decomposition on the data matrix, 
a state space triple [A, B, C] is realized to match the transfer functions 
of the system. Both discrete-time and continuous- time models are considered 
for comparisons with other methods. 

Consider a linear, time invariant system initially at rest with an m- 
dimensional input signal u(t) and an JL-dimensional output y(t), subject to 
an additive disturbance n(t), 


y(t) = 1 Y(k)u(t-k) + n(t) ; t=0,1,2,... (78) 

k=0 

where Y(k) is the impulse response matrix. This system has the transfer 
function as shown in Eq. (76) 

” -jkAtai. 

Y(j«o.) = I Y(k)e 1 At ; i-0,1,2,... (79) 

1 k=0 

The identification problem now is the following. Generate and/or 
measure an input signal u(k) (k=0,...,N) and measure the corresponding 
output signal y(k) (k=0,...,N). Based on these measurements, form an 
estimate of the transfer function, subject to the additive noise n(t), 

N -jkAtco. N -JkAtuj. 

Y ( jw. )= l Y(k)e At= l CA Be At = CG(ju> )B; i=0,...,N (80) 

0 1 k=0 ' k=0 
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N ^ -jkAtw^ 

where G(joj.) = £ he At. Our aim is that, given an estimated 

1 k=0 

transfer function Y Q (juj.), find a triple [A, B, C] of minimum order such 

that identities of Eq.(80) hold. The triple with minimum order will 
minimize the effect of noise on the identified modal parameters. 

Define a shifted transfer function by 

N -jkAtuj N -jkAto). 

Y (jo> )= l Y(k+x)e At= l CA T Be At=CA T G( jto. )Bj x-0, 1 ,2, . .(81 ) 

T 1 k=0 k=0 


A special recursive formula to compute the shifted transfer function from 
Yq ( juk ) has been developed in Ref. 26. The idea of the shifted transfer 

function is derived from the basic concept of system realization, i.e. 
Hankel matrix. 


The Eigensystem Realization Algorithm in frequency-domain begins by 
forming the rx(N+1) complex block matrix 


H (k) = 
g 


v j “o> 


V J “i> 


v j v 




Y k+t ‘jy Y k+t 0“^ ••• Y k+t 

r-1 r-1 r-1 


= V A k W 


(82) 


where 


CA 


CA 


'r-1 


, W = [G(ja) Q )B G(jw^)B ... G(jw^)B] 


and t ( i = 1 , . . . , r- 1 ) are arbitrary integers 


The matrix V is the 


observability matrix and the matrix W is the controllability matrix in 
frequency domain. Now find the singular value decomposition for the matrix 

V 0) 


* 

H (0) = P D Q ; * = complex conjugate transpose (83) 

g nnn 


where P n and Q n are orthonormal matrices in complex domain, and D n is a 

diagonal matrix with positive elements [d^ , d^, referred to as 

singular values of H (0). The rank of H (0) is determined by testing the 

S g 

singular values for zero. The pseudo-inverse of the matrix H (0) can then 

g 

be given by 

H # “ Q D _1 P* (84) 

g nnn 
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Now observe that, from Eq.(82), 

H (0) = H (O)hJh (0) = VWH^VW (85) 

o o o o o 

which implies 


WH # V = I 
g n 

where I is an identity matrix of order n 


(86) 


Define 0^ as a null matrix of order X,, 1^ as identity matrix of order p 
T 

and E^ = [0^, . . . ,1^, . . . ,0^] where 1^ is located at the ith position. With 

the aid of Eqs. (83) - (86), a minimum order-realization can be obtained 
from 




E„,H (k)E . 
£1 g mi 


T k 

= E„,VA WE , 
£1 mi 


■ E L Vt “ H g v]Ak[WH 6 VJWE mi 

' E I, H g (0)H g VAkwH g H g <0)E mi 

_T ' 1 /2 r -1 /2 * k Iir . -1/2-,_ 1 /2 # 

= E n P D. [D P VA WQ D ]D Q E . 

£1 n n n n n n n n mi 

T 1/2 -1/2 * -1/2 k 1/2 * 

= E P.d! ID * r H ( 1 ) Q D J V E , 

£1 n n n ng.nn n nmi 


(87) 


Examination of Eqs. (81) and (87) shows that the triple 


- 1/2 * - 1/2 
D Wd P H ( 1 ) Q D 
n n g n n 


-1 1/2 * 

, B = G , ( jo) . )D Q_E 


n mi 


E £1 P n D n 


1/2 


( 88 ) 


is a minimum realization derived from frequency domain analysis. To compute 

the matrix B, the integer i can be any value from 0 to N and G ^(jio..) can be 

approximated by [jukI - (lnA)/At] (see Eqs. (8) & (75)). For simplicity, it 

is set to i = 0. Based on the special characteristics of the forms B and C, 
accuracy indicators were developed in Ref. 39 to quantify the system and 
noise modes. Although the transfer function matrices Y.(ju>.) (i = 0, ..., 

. K X 

N) are in the complex domain, the block matrix H(k) can be implemented in 
the real domain by putting the real part of each individual matrix Y (ju>.) 

K 1 

in one block and its imaginary part in a consecutive block. In doing this 
way, all the computations required for the system realization become real 
arithmetic. 


Eq. (88) is developed using the discrete-time dynamic model as the 
basis so that the realized matrix A represents the state transition matrix 
(see Eq.(9)). The question arises whether a realization can be derived 

using the continuous-time model (Eq.(3)) to identify a state matrix A 
directly. Recall Eq. (76) 

N - -juj kAt 

Y (Ju.) = l Ce A Be At = CG(Jw.)B ; 1-0,..., N (89) 

U 1 k=0 


2H 



Redefine the shifted transfer functions Y^ja^) with the aid of Eq. (77) as 

Y (Ju )-(J«.) T Y (Ju.) - l Y (k) (t=0) ( jw) 1_1 _k =CA T G( joi )§; i=0 N (90) 

T 1 1 1 k=0 

The shifted transfer functions in this case are obtained by differentiating 
the impulse response function matrix t times. For example, if deflection 
sensors are used for dynamic measurements, Yq(Jio.) represent the transfer 

function matrix for the deflection responses, yjw^) describe the velocity 

responses and Y^juk) describe the acceleration responses. It is generally 

inadvisable to differentiate a measurement signal because noise effects are 
magnified. On the other hand, if accelerometers are used for the dynamic 

measurements, Y 2 (jok) represents the corresponding transfer function 

matrices. Then matrices ^ (juk) and Y 0 (juk) can be obtained by integrating 

the acceleration signals once and twice to respectively describe the 
velocity and deflection signals. The information of high frequency modes 
may be lost due to the integration process, however. 

Now substituting these shifted transfer functions (Eq. (90)) into the 
block matrix (Eq. (82) and performing the same procedures (Eqs. (83) - (88)) 
such as singular value decomposition, etc., a minimum realization identical 
to Eq. (88) will be obtained except that the state transition matrix A is 

replaced by the state matrix A. The detailed description is omitted. 
Instead, a simple example is discussed. 


Let the block matrix H (0) be formed by 

6 

Y 0^ w 0^ "■ Y 0^ W N^ 


H (0) = 
g 


^uwq) yjy ••• yjy 


= V w 


(91) 


where V = 


CA 


and W = [G(jw 0 )B GtjyB ... G( jm^)B] 


Note that the controllability matrix W is identical to the one shown in Eq. 
(82) except that B is replaced by B. Observe that 


Rgd) 


yj» 0 ) yju,) ••• yjy 

Y 2 (ju) 0 ) Y 2 (j Ul ) ... Y 2 (> n ) 


V A W 


(92) 


and assume that the number 

order of the system. Find 

* 

matrix H (0) = P D Q 
g n n n 


of rows for the matrix H (0) is greater than the 

8 

the singular value decomposition for the block 

(93) 


25 


i 



which is similar to Eq. (83). The triple 


- 1/2 #- - 1/2 -1 
A = [d’^PHIDQD 1 ^] , B = 
n n g . n n 


' 1 ( J“i )D n / 2 Q n E pi* 


T 

E ,P 
ml n n 


1/2 


(9M) 


is a direct minimum realization for the continuous-time model (Eqs.(3) & 
(4)). This simple example will be used to derive and discuss with other 
methods. 


III-3. The Polyreference Technique in Frequency Domain (Refs. 27-28) 

This subsection presents the Polyreference technique in frequency- 
domain for modal parameter identification using system realization theory. 
The methods make use of a set of transfer function matrices and shifted 
transfer function matrices to form a complex Hankel-like data matrix. By 
employing singular value decomposition of the data matrix, an orthonormal 
matrix is computed to derive an observable canonical-form realization which 
is in parallel to the canonical-form realization obtained for the 
Polyreference technique in time-domain. 


The first part of this subsection will show the direct canonical-form 
realization of the discrete-time model. The close relationship of the time- 
domain and frequency-domain Polyreference techniques is established. The 
second part of this subsection will show the -diree-t observable-form 
realization of the continuous-time model. The close relationship between 
this approach and other existing approaches (Refs. 27 & 28)— is discussed. 

Form the (r+1)x(s+1) block Hankel-like matrix 


H (0) 
g 


Y o^“o^ Y o^ w i^ *** Y o^“V 

yju> 0 ) yjy ••• Y i (j V 

yjv VJV ••• VJV 


= V w 


(95) 


where 


C = 


C 

CA 

CA r 


, W = [G(Ju 0 )B G(J«o 1 )B ... G(ju )B] 


(96) 


and r is an integer chosen such that the number £r (S. = the number of 

outputs) is greater than the order of the system. Using the same procedures 

shown in Eqs.i^) - (M8), first find the singular value decomposition of 

H (0) such that 
g 


r p t i 

H( 0) [Q Q ] = 

i 

0 

c 

Q 

1 

n 

n o 


T 

P 


0 0 _ 

o 




( 97 ) 


which yields 
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p t vwq 

o n 


»> PI = 0 
o 


( 98 ) 


T T 

where P , Q , P , and Q are orthonormal matrices and D is a diagonal 
n n o o n 

matrix. The integer n is determined by the characteristics of the system 

noise as discussed in Ref. 25. All singular values after n are considered 

as zero values. 

T T T T 

Now partition the matrices P such that P = [P P .... P ] and 

o o oO ol or 

choose square matrices P . of order i from matrices P . (i = 0, 1 

01 01 

r). Substitution of the observability matrix defined in Eq. (96) into Eq. 
(98) with the partitioned matrices P Qi produces 

r-1 „ . , . r-1 


- Z [P T ] _1 P j C A 1 = C A P = Z P T C A 1 
i=0 or 01 i=0 01 


(99) 


- T — T — 1 — T 

where P . = [P ] P . 

oi or oi 

Eqs. (99) can be rearranged into companion matrix form as 


c 


0 

X * 

0 

0 



C 

CA 


0 

0 

L l *•* 

0 



CA 

ca ; - 2 

A = 

0 

0 

0 




CA - 2 

ca ^ 1 


1 

1 

J 

O H 
O 

oi 

-P T 

02 

-P T 

^o(r- 

1) - 


r-1 

CA 

matrix 

I is an identity 

X* 

matrix of 

order 

i. 

Now it i 


( 100 ) 


Eq. (100) that the triple 


A= 


0 

£ «, 

0 

0 



0 

0 

L l 

0 



: 

: 

: 

: 

, G( jo)^ ) B= 

: 

0 

0 

0 

x x, 


•H 

3 

<\J 

1 

L, 

>H 

1 

l 

O H 

o 

-P' 1 , 

oi 

-P T 

^02 •” 

-P T 

K o(r-1 ) _ 




, C=[I £ 0 ... 0 0] 


( 101 ) 


is an ir-diraensional realization of the system. Indeed, it can be readily 
verified that 


Y 0 (j Wi ) = CG(ju).)B, Y^jw.) = CAG( joj^)B, .... Y r (ju>.) = CA^Uto^B. 


Because of the structure of A and C, it is easy to show that the 
realization is observable. However, it is not in general controllable. 
This realization is called an observable canonical-form realization. 
Therefore, it is not a minimum realization. Since the controllability matrix 
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W defined in Eq. (96) does not have the form which explicitly involves the 
state transition matrix as the observability matrix does, there does not 
exist an explicit controllable-form realization as that for the time-domain 
Polyreference technique. 

Eq. (101) is developed using the discrete-time dynamic model as the 

basis such that the realized matrix A represents the state transition 
matrix. It is natural to question whether a realization similar to Eq. 
(101) can be derived directly using the continuous-time model. The answer 
is affirmative. Replacing the shifted transfer function matrices in Eq . 
(95) by the ones shown in Eq. (90) and performing the same procedures 
(Eqs.(96) - (101)), an observable canonical- form realization identical to 

Eq. (101) will be obtained except that the triple [A, B, C] is replaced by 

[A, B, C] where A represents the state matrix for the continous-time model. 
The detailed description is omitted. Instead, a simple example is discussed 
for comparison with other existing methods (Refs. 40 - 41). 


H g (°) 


Let the Hankel-like matrix H (0) be formed as 

g 


V J “o> 

» — • 
3 
*•”> 

O 

>H 

... Y 0 (ju N ) 



O 

3 

*n 

yjy 

• • • Y ,| ( j ) 

= V w 

(102) 

o 

3 

OJ 

yjy 

• • ♦ Yg( j 0)^ ) 




where 


C 



, W = [G(j« 0 )B G(ja» 1 )B ... G(ju> N )B] 


( 103 ) 


Assume that the number of rows of the matrix H (0) is greater than the order 

& 

of the system. Following the same procedures shown from Eqs. (97) to Eq. 
(100) and using the same notations produces a polynomial equation 


-2 - T - - T 

C A + P C A + P : C - 0 
ol oO 

Eqs. (99) can be rearranged into companion matrix form as 


' C ‘ 


0 



“ C 

CA _ 

A » 

1 

1 

*Oi 
O H 
O 

-p 1 , 

ol 


. _ 


The matrix I is an identity matrix of order l. The triple 

X* 


A = 

0 


G(ja) )B = 

t 1 



l 

1 

W 
O H 
O 

ol 




(104) 


(105) 


( 106 ) 
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is a 2£-dimensional realization of the system. Examination of Eqs. (2) and 
(106) reveals that the mass matrix M, the damping matrix D and the stiffness 
matrix k for a finite-dimensional system can be related by 

M _1 K = P T . and M _1 D = P T , (107) 

oo oi 

Modal parameters of the estimated system can be obtained by solving the 
eigenvalue problem of Eq. (106). If real and undamped modes are sought, 

-T 

just solve the eigenvalue problem with the absence of the matrix P . in the 

ol 

state matrix A. The approach used here to obtain the observable canonical- 
form realization is different from that shown in Ref. 40 which first 

introduced the companion-form state matrix A shown in Eq. (106) for modal 
parameter identification.. However, both methods which use the same transfer 
function matrices to build the state space model are conceptually similar. 

The orthonormal matrix P q is computed through the application of the 

singular value decomposition to realize a companion-form state matrix. 
Since the orthonormal matrix is very close to the unity matrix, the method 
presented here thus generates a computationally well-behaved realization. 


Now consider a simplest case where the Hankel-like matrix H (0) is 

8 

formed by only two series of transfer function matrices YqUu^) and Y ^ ( j to ± ) 

(i = 0, 1, ..., N) as shown in Eq.(92). Assume that the number 2i of rows 
of the Hankel-like matrix is at least twice the order of the system. 
Following the same procedures as shown from Eqs. (97) to (99) yields a first 
degree polynomial equation as 

CA + pIc=0 : (108) 

oO 

The triple 

A = -PJ Q , G(j(u.)B = Y 0 (j W .) , C - I t (109) 


is a ^-dimensional realization of the system. Eq. (108) was introduced in 
Ref. 41 for modal parameter identification for flexible structures. In 
contrast to the minimum realization shown in Eq. (94), this realization (Eq. 
(109)) is not of minimun order in the sense that the identified state matrix 
is usually oversized if the order of the system is not known a priori. 
Furthermore, the number of sensors L must be greater than the order of the 
system. 


III-4. An Alternate method for the ERA-FD and the Polyreference Technique in 
Frequency Domain 

A minimum order of canonical-form realization in frequency domain is 
generally impossible for multi-input and multi-output systems due to the 
constraint that the realized state matrix is a companion form. If the 
constraint is released, a minimum order realization in frequency-domain can 
be obtained from Eq . (101). The procedures to derive the minimum order 
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realization are identical to that shown in Section III-5 for the time domain 
case. Using the block Hankel-like matrix (Eq.(95)) and notations defined in 
Section 11-5, the triple 


D n 1/2CE ir P n 3#C °iL E !Lr P n 3D i /2 ' 
n Hr n x, x,r n n 


-1 1/2 * 

B - G ( ju> . )D Q E , 
in n m 


T 1 /2 

C = Wn (110) 


is a minimum realization of order n. Here, [o.E. P ] simply means the 

X Xj l 11 

matrix obtained by deleting the last fc rows of the matrix P n and [E^ P n ] 

represents the matrix obtained by deleting the first l rows of the matrix 

P . Since P is an orthonormal matrix, a special and efficient procedure 
n n 

can be developed to compute the pseudoinverse of the matrix E^P using the 
matrix inversion lemma (Ref. 19). 

Based on Eq. (108), a close link between between the ERA-FD and the 
Polyrefernce technique in Frequency domain is established. A minimum 
realization similar to Eq. (108) can also be derived for the direct 
realization of the continuous-time model. 


After several methods in f requency- doma i n for modal parameter 
identification are derived using system realization theroy, the derivation 
of the least squares regression technique in frequency domain becomes 
trivial and thus is omitted. Following the same procedure as shown in the 
section II- 6, the reader is encouraged to derive it on his own. 
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IV. CONCLUDING REMARKS 


In this report, several methods for modal parameter identification have 
been presented and derived using system realization theory. The relations 
between different techniques are reasonably well understood and the choice 
of methods can be done largely on the basis of the final purpose of the 
identification, for example control of flexible structures. Most methods 
are claimed to work well on simulated and test data. In spite of a large 
body of literature on identification, there are few papers which compare 
different techniques using experimental data. Unfortunately, conclusive 
results have not been obtained. It is more fun to dream up new methods 
than to work with somebody else's scheme. However, for a person engaged in 
application, it would be highly desirable to have comparisons available. 
This report illustrates the mathematical relations among several recently 
developed methods via system realization theory, which provides a basis and 
insight for comparison and evaluation. The practicing engineer would still 
not be satisfied, however. It would be necessary to have a selection of 
test data from real systems to which several techniques are tried and 
compared for evaluation. 

It is hoped and expected, through combined efforts in the control and 
structural dynamics disciplines, that the field of modal parameter 
identification is moving towards more unification and that there will be 
more comparisions of different methods. Two purposes of this report is to 
contribute to the goal of unification and to serve as a starting point to 
stimulate more research toward this goal. It is believed that interaction 
with other disciplines such as controls, artificial intelligence, etc. is 
essential for accomplishment. 

For a report like this size, it is difficult to be complete. There are 
many other techniques available in the fields of modal parameter 
identification. The reader is directed to the References, Bibliography on 
System Realization Theory, and other literature for further information. 
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